Speeding up disease extinction with a limited amount of vaccine 
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We consider optimal vaccination protocol where the vaccine is in short supply. In this case, disease 
extinction results from a large and rare fluctuation. We show that the probability of such fluctuation 
can be exponentially increased by vaccination. For periodic vaccination with fixed average rate, the 
optimal vaccination protocol is model independent and presents a sequence of short pulses. The 
effect of vaccination can be resonantly enhanced if the pulse period coincides with the characteristic 
period of the disease dynamics or its multiples. This resonant effect is illustrated using a simple 
epidemic model. If the system is periodically modulated, the pulses must be synchronized with the 
modulation, whereas in the case of a wrong phase the vaccination can lead to a negative result. 
The analysis is based on the theory of fluctuation-induced population extinction in periodically 
modulated systems that we develop. 
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I. INTRODUCTION 

Spreading of an infectious disease is a random pro- 
cess. An important source of the randomness is demo- 
graphic noise associated with the stochastic character of 
the events of infection, recovery, birth, death, etc. In 
a large population demographic noise is small on aver- 
age, and the infection spread leads to an endemic state 
where a certain fraction of the population stays infected 
for a long time. The disease, however, can disappear as 
a result of a large rare fluctuation, an unlikely chain of 
events where, for example, susceptible individuals hap- 
pen to avoid getting infected while infected ones recover 
[l|-l3|. Then, if there is no influx of infected individu- 
als from the outside, the population will be disease-free. 
Such spontaneous disappearance of a disease is an exam- 
ple of population extinction studied in stochastic popu- 
lation dynamics. 

Spontaneous extinction is also important for various 
physical and chemical reaction systems. This is a conse- 
quence of the underlying similarity of the dynamics that 
result from short random events, like collisions between 
molecules that lead to chemical reactions and interac- 
tions between individuals that lead to the disease spread. 
Extinction can be understood for different types of sys- 
tems within the same general formalism, which provides 
a broader scope for the present paper. Moreover, the 
method of optimal control of extinction that we propose 
can be applied to systems of various types. 

A conventional way of fighting epidemics is via vac- 
cination. If there is enough vaccine, the infection can 
be eradicated "deterministically" by eliminating the en- 
demic state [1]. The amount of available vaccine, how- 
ever, is often insufficient. The vaccine may be expensive, 
or it may be dangerous to store in large amounts, as in 
the case of anthrax, or it may be effectively short-lived 
due to mutations of the infection agent, as for HIV [sl 
and influenza 6]. 

Even where the endemic state may not be eliminated 
deterministically, vaccination can dramatically affect the 



stochastic dynamics of the epidemics. The underlying 
mechanism is the change of the rate of large fluctuations 
leading to disease extinction. For a well mixed popu- 
lation, this rate We is usually exponentially small for 
a large total population size N ^ 1, We oc exp(— Q) 
with Q (X N, % 0-[il- We caU Q the disease extinc- 
tion barrier. Vaccination changes the value of Q/N. In 
turn, this changes the disease extinction rate exponen- 
tially strongly. This effect was previously discussed for 
vaccination applied at random [l3|. 

The goal of this paper is to find an optimal way of 
administering a limited amount of vaccine which would 
maximally increase the disease extinction rate. We find 
a vaccination protocol that applies for a broad class of 
epidemic models. Our approach is based on the obser- 
vation that, in a large fluctuation that leads to disease 
extinction, the population is most likely to evolve in a 
well-defined way. It moves along the most probable path 
in the space of the dynamical variables which characterize 
different sub-populations [Tsl-fisj. Vaccination perturbs 
the system as it moves along the optimal path. One can 
think of vaccination as "force" and its effect as "work" 
done on the system. This work reduces the barrier Q. 
The problem then is to maximize the work for given con- 
straints on the vaccine. 

Optimization of the effect of vaccination resembles an- 
other problem of optimal control of random systems: 
controlling large fluctuations in noise-driven dynamical 
systems by applying an external field with a given av- 
erage intensity jl6l Il7l|. There are, however, important 
differences, which come from the very nature of the con- 
trol field. Indeed, vaccination only reduces the number 
of susceptible individuals. In other words, as a control 
field, vaccination never changes sign. Then, as we find, if 
the available amount of vaccine is constrained by a given 
mean vaccination rate, the optimal vaccination protocol 
turns out to be model- independent. This applies also to 
using a limited amount of medications and other situa- 
tions where the control field drives the system only in 
one direction. 
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We assume that vaccination is apphed periodically in 
time. In this case, the optimal protocol is to apply the 
vaccine as a sequence of 5-like pulses. The disease ex- 
tinction rate can strongly depend on the period of this 
sequence. Furthermore, the extinction rate can display 
exponentially sharp peaks when the pulse sequence pe- 
riod is close to the characteristic period of oscillations 
of the system in the absence of fluctuations, or to its 
multiples. We illustrate this resonant phenomenon for 
the Susceptible- Vaccinated-Infected-Recovered (SVIR) 
model. 

Epidemics often display seasonal modulation [l8|. It 
is natural to apply a vaccine with period equal to the 
modulation period. As we show, there is a qualitative 
difference between the effect of a periodic vaccination in 
this case and in the case where seasonal modulation is 
absent. For a system with seasonal modulation, an im- 
properly applied pulsed vaccination can actually reduce 
the disease extinction rate and therefore prolong the du- 
ration of the epidemic. The overall effect of the pulsed 
vaccination critically depends here on the phase at which 
the periodic pulses are applied. 

The analysis of periodic vaccination, with and with- 
out seasonal variations, necessitates a general formula- 
tion of the extinction problem in periodically modulated 
stochastic populations. We extend the previous results 
for single-population systems [l^ [l^l to provide a com- 
plete extinction theory for multi-population systems in 
the eikonal approximation, and emphasize the distinction 
from the well-understood problem of switching between 
metastable states in periodically modulated systems with 
noise ^)\. 

Section II describes the class of epidemic models we 
consider in this work and develops an eikonal theory 
of disease extinction rate in periodically modulated sys- 
tems. Section III formulates the optimization problem for 
vaccination and presents its solution for a time-periodic 
vaccination in the limit of a small average vaccination 
rate. In Section IV we discuss the vaccination-induced 
reduction of the disease extinction barrier for two types 
of constraints on the vaccination period, a limited life- 
time of the vaccine and a limited vaccine accumulation. 
Section V illustrates, on the example of the stochastic 
SVIR model, the phenomenon of resonant response to 
vaccination. Section V contains concluding remarks. 



II. THE DISEASE EXTINCTION RATE 

A. The model of population dynamics 

We consider stochastic disease dynamics in a well- 
mixed population which includes infected (/) and sus- 
ceptible (5) individuals and possibly other population 
groups such as recovered or vaccinated. The system 
state is described by a vector X — {S,I, . . .) with in- 
teger components equal to the sizes of different popula- 
tion groups. Along with X it is convenient to consider a 



quasi-continuous vector x = X/A^, where N is the char- 
acteristic total population size, ^ 1. We assume that 
the population dynamics is Markovian. It is quite gener- 
ally described by the master equation for the probability 
distribution P(X,i), 

P{X,t) = ^[M/(X-r,r,t)P(X-r,t) 

r 

- W{X,r,t)P{X,t)]. (1) 

Here W(X.,r,t) is the rate of an elementary transition 
X X + r in which the population size changes by r = 
(ri,r2,...). Examples of such transitions are infection 
of a susceptible individual as a result of contacting an 
infected individual, recovery of an infected individual or 
arrival of a susceptible individual. 

We assume that there is no influx of infected individ- 
uals into the population. Therefore, there are no tran- 
sitions from states where there are no infected to states 
where infected are present, 

VK(X,r,t)=0 for = 0, ^ 0, (2) 

where subscript E is used for the component of X which 
enumerates infected, Xe = I- 

In the neglect of demographic noise the population dy- 
namics can be described by the deterministic (mean-field) 
equation for the population size x, 

i = ^rw;(x,r,0, w{^,v,t) = W{X,v,t)/N. (3) 

r 

It immediately follows from Eq. ([l} if the width of the 
probability distribution P(X, t) is set equal to zero. 

1. Stationary systems 

We start with the case where the transition rates 
W{X, r) are independent of time. An endemic state, 
where a finite fraction of population is infected for a long 
time, corresponds to an attracting fixed point x^i of the 
dynamical system, Eq. ([3]). We will assume throughout 
this work that there is only one such point. We will also 
assume that Eq. (j3|) has one fixed point xg in the hy- 
perplane xe — The state is stable with respect to 
all variables except Xe- We call it the disease extinction 
state. If Xe > Q (there is a nonzero number of infected), 
the deterministic trajectory leaves the vicinity of and 
approaches the endemic state xa. 

Due to demographic noise the endemic state is ac- 
tually metastable. A rare large fiuctuation ultimately 
drives the population into a disease-free state. The most 
probable fluctuation brings the system to the fixed point 

El, 

IT^ . The probability of such a fluctuation per 
unit time, i.e., the disease extinction rate We, is given 
by the probability current to xg, similarly to the prob- 
lem of escape from a metastable state 22]. For time- 
independent WiX, r) this current is quasistationary for 
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times tr '^t <^ W~\ where tr is the characteristic relax- 
ation time for the noise- free motion described by Eq. ^ . 

We note that, even though the state xg is a saddle 
point in the mean-field approximation, it differs from 
the saddle-point states encountered in the problem of 
switching between metastable states of reaction systems. 
In the case of interstate switching, the rates of elemen- 
tary transitions W(X, r) in the unstable direction are 
nonzero, and ultimately fluctuations drive the system 
away from the saddle point. In the case of extinction, 
fluctuations around X5 occur only in the extinction hy- 
perplane, whereas the probability of exiting this hyper- 
plane is zero. 

If the system has, in the mean-field approximation, 
more than one steady state away from the extinction hy- 
perplane, extinction can go in steps: from the endemic 
state to another steady state and, ultimately to the ex- 
tinction hyperplane. In particular, if the only additional 
steady state is a saddle point at the boundary between 
the basins of attraction of x^i and xg, the problem of 
extinction can be reduced to the problem of escape over 
this saddle point [23l |. 

2. Periodically modulated systems 

The above picture can be extended to the case where 
the transition rates are periodic functions of time, 
W{X,r,t + T) = W(X,r,t). Periodicity of some of 
the rates in time is a natural way of modeling sea- 
sonal variations of epidemics The attracting so- 
lution of Eqs. which describes the endemic state, 
is no longer stationary. We will assume that this solu- 
tion, XA{t), is periodic in time with the same period T, 
x^(t + T) = x^(i). The asymptotic disease extinction 
state xs{t) is also periodic in time; it lies in the hyper- 
plane xe — 0- 

The most probable fluctuation which causes extinction 
of the disease corresponds to a transition from XA(i) to 
xs{t) [2^. An important characteristic of such a tran- 
sition is the period-averaged disease extinction rate We- 
lt can be introduced if the modulation period T <^ 
and, in addition, tr ^ . In this case, for time t such 
that tr,T <$:_ t <^ a quasi-stationary time-periodic 

probability distribution is formed, centered at XA{t). The 
probability current from x^ {t) to x^ (t) is also periodic in 
time, and the period-averaged value of this current gives 
We [20| , in a direct analogy with the problem of switch- 
ing between metastable states in noise-driven dynamical 
systems (H-Ill. 

B. Eikonal approximation 

1. Equations of motion 

We will be interested in evaluating the disease extinc- 
tion barrier Q which gives the exponent in the disease 



extinction rate. We oc exp(— Q), at N ^ 1. This bar- 
rier is entropic in nature, as it results from an unlikely 
sequence of elementary transitions. It can be found by 
either solving the mean first passage time problem for 
reaching X5 (t) Jj-tB^ or by calculating the tail of the quasi- 
stationary probability distribution P(X, t) for x close to 
xs(i) jlll. [13I . Here we choose the latter strategy and 
determine, to the leading order in N, the logarithm of 
the distribution tail. We seek the solution of Eq. ^ in 
eikonal form, P{X,t) = exp[-Ns{x,t)] (H-iOl. In the 
limit of large N, from Eq. ^ we obtain the following 
equation for s(x, i): 

dts = -H{^,d^s,t), (4) 
iJ(x, p, = V w(x, r, t) [exp(pr) - 1] . 

^ — 

Here, we have taken into account that, typically, |r| ^ N, 
and VF(X, r, t) depends on X polynomially, whereas P is 
exponential in X. Therefore we expanded P(X + r,t) w 
P(X, t) exp{—rdxiS) and replaced, to the leading order in 
u;(x - N-^r,r,t) by w{x,r,t). 
Equation ^ has the form of the Hamilton- Jacobi 
equation for an auxiliary Hamiltonian system with 
Hamiltonian 7?(x, p,i); s(x, i) is the action of this sys- 
tem. The Hamilton equations of motion are 

X = V rw(x,r,i)eP"", 

z — rf-r 

P = - V ax^i;(x,r,i)(eP'--l). (5) 

— 

These trajectories determine, in turn, the most probable, 
or optimal, trajectories that the system follows in a fluc- 
tuation to a given state x at time t. We will calculate 
action s{x,t) using these trajectories and thus find the 
exponent in the distribution P(X,i). 

2. Boundary conditions for the optimal extinction 
trajectory 

To find the boundary conditions for Hamiltonian tra- 
jectories ([U, we note that the quasi-stationary distribu- 
tion P(X,t) has a Gaussian maximum at XA(i). This 
means that, close to attractor x^(t), action s(x, i) is 
quadratic in |x — xa)\ for stationary systems, whereas 
for periodically modulated systems s{x,t) = s(x, t -I- T) 
is quadratic in the distance from trajectory x^(t) 31]. 
On the Hamiltonian trajectories that give such action, 
the momentum p = — >■ for x — J> x^(t), and since 
X = X.A (t),p — is a fixed point (a periodic trajectory) 
of the Hamiltonian dynamics, the trajectories of interest 
start at t ^ —00, 

s(x,t) = / dtL{±,x,t), (6) 

J —00 

L(x, X, i) - V z«(x, r, t) [(pr - l)eP'- + 1] . 

^ — 

In the Lagrangian L, Eq. ([6]), p should be expressed in 
terms of x, x using Eq. ([S|). Since w > 0, we have L > 0. 
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Extinction barrier Q is determined by Ns{yi,t) for x 
in the extinction hyperplane, xe = 0- In the spirit of the 
eikonal approximation, we have to find such (x, t) in this 
hyperplane that s(x, i) be minimal. The minimum de- 
termines the boundary conditions for the optimal Hamil- 
tonian trajectory of extinction, (xopt(i), Popt(^))- The 
condition that s(x, is minimal with respect to Xi^E 
on the extinction hyperplane means that pi = dx^s — 
ioT i ^ E as the trajectory (xopt(t), Popt(i)) approaches 
the hyperplane. The minimum of s(x, t) with respect to t 
within the period of modulation is reached if i?(x, p, t) 
as the trajectory approaches the hyperplane. 

A consequence of conditions iJ(x, p,t) — i> and 
Pi^E — > is that the momentum component pE remains 
bounded on trajectory (xopt(i), Popt(i))- Indeed, near 
the extinction hyperplane, ^ 1, we have 

ifi = r£;w(x, r, t)eP'" 

L — 

-xeY] rE[dw{^,v,t)/dxE]x^^oe^^'^^. (7) 

Here, we assumed that w(x, r, t) is nonsingular at xe — > 
and, since w = for a;^ = and r^; ^ [cf. Eq. ([5])], 
we expanded w iw xe to the lowest order. Let us assume 
now that \pe\ ^ oo for x^; — > 0. Then only the term with 
maximal —rE = —rEm should be kept in the sum over 
r^; in Eq. ([7]); it is also clear that pe should be negative, 
otherwise the trajectory would not approach xe = 0. 
From the Hamilton equation for pE and Ec^. (O it follows 
that dpE/dxE ~ XEVEm- This relation, along with 
the explicit form of the Hamiltonian H, show that, '\ipE 
were diverging for a;^; — >■ 0, the Hamiltonian would not 
become equal to zero but would remain w dw/dxE with 
the derivative calculated for xe = ^ and = rEm- This 
contradiction shows that the assumption \pe\ — > oo is 
wrong, Pe remains limited for xe ^ Q- 

Equation^ shows that xe — > exponentially as 
t — > oo. As a;^; approaches zero, variables Xi^E are 
approaching the equilibrium position in the hyperplane 
Xe = ^- This happens because Pi^E and the dy- 
namics of Xi^E in the hyperplane is described by the 
mean- field equations, Eq. ([3]). Therefore, 

/oo 
dtL{±,jc,t), (8) 
-oo 

x{t) xs{t), p{t) Ps{t) for t ^ oo. 

Function Ps{t) is periodic in time, with Pi^E = and 
with hitherto unknown component pE{t), which is dis- 
cussed below. 

The optimal trajectory (Xopt(i), Popt(i)) is a hetero- 
clinic Hamiltonian trajectory that goes from periodic or- 
bit (xyi(t),p = 0) to periodic orbit {xs{t),ps), and the 
action for extinction s^xt is calculated along this trajec- 
tory. The trajectory Xopt(i) is the optimal path to disease 
extinction: it describes the most probable sequence of el- 
ementary transitions leading to extinction. We note that, 
in periodically modulated systems, there is one optimal 
path per period, whereas in stationary systems trajecto- 
ries (xopt(i), Popt(O) are time-translation invariant. 



3. The t ^ oo value of the momentum on the Hamilton 
trajectory 

The momentum component pe is generically nonzero, 
as found for stationary systems [ij, [s^l ■ For peri- 
odically modulated systems, one can show that »e 7^ 
by extending the arguments presented in Ref. [IJ, |33 |. 
This amounts to showing that the stable manifold of the 
periodic orbit (x5(t), p = 0) lies entirely in the invariant 
hyperplane a;^; = 0, Pi^E — 0, and, as a consequence, 
does not intersect the unstable manifold of the periodic 
orbit {x.A{t),p = 0) . Such intersection is necessary in 
order to have a heteroclinic trajectory that would go from 
(xA(t),p = 0) to (X5W,P = 0). 

The hyperplane xe = 0, Pi^E = is formed by trajec- 
tories 

x^-^E = ^Jw(x,r,t)]^^^orj, (9) 
PE = - J2j9.M^,r,t)l,^,{eP-^- -I). 

The invariance of this hyperplane is a consequence of 
Eq. which leads to Pi^E — and xe = for Pi^E — 
and a;^; = 0. 

To prove that the stable manifold of {xs{t), p = 0) lies 
entirely in the invariant hyperplane xe = 0, Pi^E = 0, we 
first show that the trajectories, which are described by 
Eq. ^ and which start close to the state (x5(i), p = 0), 
approach this state for t — > 00. Then, since the dimen- 
sion of the hyperplane xe — 0, Pi^E = is equal to 
the dimension of the stable manifold of {xs{t),p — 0), 
we conclude that the stable manifold indeed lies in the 
hyperplane. 

Equations ^ for Xi^E sue the mean- field equations 
in the extinction hyperplane xe = 0, cf. Eqs. and 
therefore Xi — )■ (x5(t))i for t — > 00. Linearization of 
Eq. ^ tor Pe about {xs{t),p — 0) gives 

PE = -J2j9^Bwi:ii,r,t)]^^f^^^pErE- (10) 

We compare this equation with the mean-field equa- 
tion for Xe near :x.s{t). The latter has the form xe = 
XEj2r'^E[dxEw{:>c,r,t)]-^^(^ty Since the state X5(t) is 
unstable in a;^; direction in the mean- field approxima- 
tion, from Eq. pO)) pe/pe < 0. Therefore, all trajecto- 
ries on the hyperplane xe — ^, Pi^E = close to the 
state (x5 (t) , p — 0) approach this state asymptotically 
ast ^ 00, and thus the stable manifold of ((x.s{t),p ~ 0) 
lies in this hyperplane. 

From the above analysis one concludes that there 
are no Hamiltonian trajectories that would go from 
(xA(i),p = 0) to (x5(t),p = 0). Therefore the opti- 
mal trajectory leading to extinction should go to a state 

(xs(t), [ps{t)]E) with [ps{t)]E + 0. 
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III. OPTIMAL VACCINATION 



as 



A. Constraint on the vaccination protocol 

Vaccination increases the number of individuals who 
are at least temporarily immune to the disease. It thus 
reduces the pool of susceptible individuals and ultimately 
leads to a reduction of the number of infected. When the 
available amount of vaccine is small, so that the disease 
extinction still requires a large fluctuation, the goal of 
vaccination is to reduce the disease extinction barrier Q. 

An outcome of vaccination is often modeled as the cre- 
ation of a sub-population of vaccinated individuals out 
of susceptibles. The corresponding elementary transi- 
tion rate is W^(X,r) = S,Q{t)Xs for rs = —1, ry = 1 
and ri^sy = Oj where subscripts V and S refer to 
vaccinated and susceptible individuals, respectively, and 
^o(i) is the control field that characterizes the vaccina- 
tion (subscript S should not be confused with subscript 
S used to indicate the extinction state). Another model 
is vaccination of newly arrived susceptibles, which leads 
to an effective reduction of the arrival rate ijlN . In this 
model, the elementary transition rate for the arrival is 
I^(X,r) = N[^l - io{t)] for rg = 1 and n^s = 0, with 
£,o{t)N being the change in the arrival rate due to vacci- 
nation. 

We will consider a general model where vaccination 
modifies the rate of an elementary transition of a certain 
type; the change of the population in the corresponding 
transition is r^. The field £,o{t) characterizing the vacci- 
nation is assumed to be weak. The affected rate has the 
form W{X,r^,t) = W'^^l {X,r^,t) + (oit)W^^\X,r^,t), 
with II/^") being the rate without vaccination. The vac- 
cination either increases or decreases the rate, as for tran- 
sitions from susceptibles to vaccinated or for vaccination 
of newly arrived susceptibles, respectively. Therefore, we 
will assume without loss of generality that £,o{t) > and 
that W^^^^(X, r^,t) is either positive or negative. We con- 
sider models in which the number of susceptibles changes 
by 1 in an elementary transition associated with vacci- 
nation, (r^)s — ±1. We note that the analysis can be 
immediately extended to describe other processes, like 
modification of the infection rates Q or recovery accel- 
eration by administrating medicine. 

It should be noted that the vaccination model adopted 
in this work is probabilistic by nature. An alternative 
is where vaccination is done in a pre-determined fashion, 
when a certain number of individuals are vaccinated per 
unit time at a given time. The analysis of such determin- 
istic vaccination lies beyond the scope of this paper. 

We will assume that vaccination is periodic, £,o{t) = 
f {t ~^ T), and that the amount of vaccine available per 
period T is limited. We model this limitation as a con- 
straint on the ensemble-averaged number of individuals 
vaccinated per period T. The constraint can be written 
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Here, S is the average vaccination rate rescaled by the 
characteristic population size N. The constraint is well- 
defined for tr <^ t W-\ where P{X,r,t T) w 
P(X,r,t). Since for ^ 1 the population distribution 
sharply peaks at the endemic state X^(i), the sum over 
X in Eq. ^ can be replaced by \w'-^\XA{t),r^,t)\, in 
the leading order in 1/N. 

In the presence of vaccination, one can still seek a so- 
lution of the master equation in the eikonal form. The 
exponent Q in the extinction rate is again given by the 
action of an auxiliary Hamiltonian system, Eq. ([5]). The 
Hamiltonian now has the form 

£r(x, p) = (x, p) + ^0 Wi?(^nx, P), 
if(°)(x,p)=^ «;(«)(x,r,t)(eP"--l), (12) 

i7(i)(x,p) = u;(i)(x,r4,i)(eP''« - 1). 

Our goal is to find the optimal form of S,o{t) which would 
minimize the disease extinction barrier Q subject to con- 
straint ([TT|) . Since w^^^ (x^(t), r^, t) is a known periodic 
function of time, we can equivalently search for the op- 
timal vaccination rate ^{t) = S,o{t) | w'^"'^^ (x^(t), r^, t) | . It 
minimizes the functional 

Jo 

at) = at + T) - Ut) \w'^'\^A{t),r^,t)\ > 0, 

where A is the Lagrange multiplier. The functional Soxt [C] 
is given by Eq. ([5]), where the Lagrangian corresponds 
to the Hamiltonian ()12p and depends on the vaccination 
rate £,{t). 

B. Vaccination protocol for stationary systems 

1. Double optimization problem 

For a low average vaccination rate S it suffices to keep 
in the action s^xt only the leading-order term in ^(t). 
Since Hamiltonian (jl2p is linear in ^, this term is linear 
in ^, too. In the spirit of the standard perturbation the- 
ory for Hamiltonian systems j33;j, the change in the ac- 
tion, caused by the small perturbation, can be calculated 
along the unperturbed trajectory {:x.l^p^{t),pl^J^{t)) of the 

Hamiltonian which describes the optimal path of 

disease extinction in the absence of vaccination. 

We start with the case of systems, which are stationary 
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in the absence of vaccination. For such systems 

«ext[eW]=4xl + 4xlKW], (14) 



s(;,l[C(t)]=min 
to 



X{t) = -if(i)(x(°) (0,pi^pUO) w^'H^A,r^ 



dixit ~to)m, 



(0) 



The quantity x(t) is called logarithmic susceptibility [ij, 
[20l . [34 [35} ; it gives the change of the logarithm of the 
extinction rate, which is linear in the vaccination rate for 
moderately low vaccination rate. 

The minimization over to in Eq- (|14p accounts for lift- 
ing the time-translational invariance of the optimal ex- 
tinction paths mentioned earlier. For £^{t) = 0, extinc- 
tion can occur at any time (i,. ^ i ^ VF~^) with rate 
We- Periodic vaccination synchronizes extinction events; 
it periodically modulates the extinction rate, and the 
modulation is exponentially strong for iVjs^xtl ^ 1 (see 
below). Formally, in a modulated system there is only 
one optimal extinction path per period, as explained in 
Sec. II, which is here reflected in the minimization over 
to. This optimal path minimizes the disease extinction 
barrier Q = A^Sext [H, S H Hi]. Equation ^ is 
closely related to the Mel'nikov theorem for dynamical 
systems [13, HI]. 

The constraint for minimizing the action over ^(<) in 
Eq. (fT3|) has a form of an integral over the vaccination 

period T. It is therefore convenient to write action Sg^t 
also in the form of such an integral, 

pT 

4ilK(t)] - mill / dtat)xT{t^to), 



^0 Jo 



XT(t)= J2 xit + nT). 



(15) 



The function xrit) is obtained by superimposing the 
parts of xit) which differ by T. By construction, XT{t) 
is periodic in time t. 



2. Temporal shape of optimal vaccination 

To find the optimal shape of vaccination rate ^(t) we 
first minimize the time integral in the variational prob- 
lem Eqs. (fT3)) - with respect to £,{t) for a given tg. 
Since ^(t) > 0, it is convenient to perform the minimiza- 
tion with respect to $.^^'^{t) rather than ^(t). The min- 
imization shows that C^^^(i) 7^ only for t = t\, where 
t\ is given by equation XT{t\ ~ to) — ~X/T. From the 
constraint on the period-averaged ^(i) we then have 



at)^ETj26{t-tx + nT). 



(16) 



Substituting this expression into the functional Soxt and 
minimizing with respect to to, we obtain the action in a 



simple explicit form 



-^cxt 



s^^{ = min xrit)- 

0<t<T 



(17) 



Alternatively, this expression can be rewritten in terms 
of the Fourier transform of the logarithmic susceptibility: 



sixt = S min V xin^) exp[innt], 

t ' * n 



(18) 



x(w) 



dtx{t) exp(ia;i). 



where = 2tt/T is the cyclic frequency of vaccination. 

We are interested in the solution for which Sg^t is nega- 
tive, which requires minxTit) < 0. Only in this case will 
vaccination reduce the barrier for disease extinction and 
increase the disease extinction rate. The barrier reduc- 
tion due to vaccination, Q^^^ = Ns'^^ oc TVS, becomes 
large for 1 even if the average vaccination rate S is 

small. Thus, for not too small vaccination rates, where 
the eikonal approximation is valid [2^ [sj] , the effect of 
vaccination on the disease extinction rate is exponentially 
strong. 

The expression for the action change Sg^t, Eq. (flT)) . can 
be also obtained in a more intuitive way. Indeed, since 
^(t) is non-negative, it follows from Eq. (fT5|) that 



■^cxt 



m] > ^}^XT{t) I dtm = SrminxT(<). (19) 
* JO ' 

The minimum is provided hy ^{t) ='E.T J2n S(t — tmin + 
nT). Formally, tmin is the instance of time where XT(t) 
is minimal. In fact, it is the optimal path that adjusts to 
the vaccination pulses so as to increase the probability of 
disease extinction. This provides the mechanism of syn- 
chronization by vaccination. Equation (|19p immediately 
leads to Eqs. and ([T7)) with t\ replaced by tmin- 

In addition to the constraint on the average vacci- 
nation rate, there may be an upper limit on the in- 
stantaneous vaccination rate, which is imposed by con- 
dition w(x, 

the case where w'^-' (x, r^, i) < 0, as for vaccination of 
newly arrived susceptibles, this condition is met provided 
Ut) < Com = min{u.(o)(x,r5)/|u;(i)(x,r5)|}. In this 
case the optimal vaccination protocol changes. 

The new protocol can be found from the variational 
problem (fT3|) by changing from ^(t) =^o(i) l^'^^^H^^; i"?) | 
to an auxiliary function r](t) such that S,o{t) = Com[l -|- 
7y^(i)]~^, and then finding the minimum of Sext with 
respect to 77 (i). This choice satisfies the constraints 
< £,o{t) < Com- Variation with respect to j]{t) shows 
that Soxt has an extremum for ri(t) = or ri(t) — 00 for 
t ^ t\, where tx is given by equation XT{t\ — to) ~ —\/T. 
The value of r](t) at the isolated instances t — t\\s arbi- 
trary. Only regions where rj{t) — 0, so that £,o(t) = ^om, 
contribute to Scxt- Obviously, Sext is minimal when 
io{t) = Com for \t - tmin\ < At/2, where tmin is the 
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time when xt (t) is minimal and At is determined by the 
average vaccination rate S. In other words, the vaccina- 
tion rate ^ {t) has the form of periodic rectangular pulses 
of width At, centered at tmin + n-T, n = 0, ±1, ±2, .... 
The pulse width is 

At= , I (f.y ^. (20) 

Since the vaccination rate is limited by the rate 
of elementary transitions without vaccine, we have 
Com |w'-^-'(x^, r)| < t^^. Then for weak vaccination, 
ST < 1, from Eq. ^ At <s: U. Therefore, xrit) = 
XT(^min) during the pulse of C(t), to the leading order in 
ST [we note that XTit) may vary on a time scale shorter 
than tr, see below; however, this time scale is always long 
compared to At for sufficiently weak modulation]. The 
resulting change of the action is again given by Eq. (fT7| . 

C. Vaccination protocol for periodically modulated 
systems 

Optimal vaccination in periodically modulated systems 
requires a separate consideration. Here, there is only one 
optimal extinction path per period T in the absence of 
vaccination. When the average vaccination rate S is low, 
vaccination with the same period T will only weakly per- 
turb this path. To first order in S, the linear in C term in 
the action still has the form of Eq. but without min- 
imization over to. Since ^(t) > 0, the minimum of action 
is still achieved for C(t) = 2T^^ S{t — tmin + nT), but 
now tniin is uniquely determined by the strong modula- 
tion. In other words, the modulation uniquely determines 
the phase of the optimal vaccination pulses. The result- 
ing expression for s^^{ for the optimal vaccination pro- 
tocol has the form of Eq. p7)) . If the vaccination pulses 
are applied at a wrong time, i.e. if the phase difference 
between the vaccination and the modulation differs from 
the optimal one, the vaccination will be not as efficient 
and may even be harmful: it may prolong the lifetime 
of the endemic state by increasing the disease extinction 
barrier Q. 



1. Limited lifetime of the vaccine 

The vaccination period T is naturally limited by the 
effective lifetime Ty of the vaccine. This lifetime is usually 
determined by the maximum storage time of the vaccine 
and/or by the mutation rates of the infectious agent. If 
Tv is long, Tv S> tr, vaccination can be made most ef- 
ficient by increasing the vaccination period up to t^. 
Indeed, as it follows from Eq. ([TT)) . Sg^t C)C T in this case. 
This result is easy to understand. Even though a de- 
crease of the vaccine pulse frequency = 2'k/T causes 
a decrease of the prefactor in the disease extinction rate 

We, the exponential factor exp(— TVSgxt) in We increases 
sharply. Indeed, it can be seen from Eqs. and 
that, as the system moves along the optimal path 
to extinction, x(t) is significant when the system is far 
from the stationary states Xyi and x^ . The characteristic 
time scale of this motion is ^ tr. For T ^ tr we have 
mmo<t<T Xrit) ~ miuf x{t) and 

si;,! = STmmx(t), t„ > T » t,. (21) 

In the opposite limit of <C t,., and thus T <^ tr, we 
have from Eq. (fT8)) 

4;,l = Sx(0), t, »r,>T, (22) 

In this case the vaccination-induced reduction of the ex- 
tinction barrier is independent of the vaccination period 
and is determined by the zero-frequency component of 
the logarithmic susceptibility. 

An interesting situation may occur in the intermediate 
range tr if, in the mean-field description, the system 
approaches the endemic state in an oscillatory manner. 
In this case function x(t) is also expected to oscillate. 
The oscillations are well-pronounced if their typical fre- 
quency is uiQ ^ tr^ . It is clear from Eq. (|18p that a strong 
effect on disease extinction can be achieved by tuning the 
vaccination frequency ft — 27r/T or its overtones in res- 
onance with ujo ■ An example of such a resonance will be 
discussed in Sec. V. 



IV. DISEASE EXTINCTION BARRIER AS A 
FUNCTION OF VACCINATION PERIOD 

The vaccination-induced reduction of the disease ex- 
tinction barrier Q^^^ = A'^s^^j, as given by Eqs. ([IT]), de- 
pends on the interrelation between the vaccination period 
T and the characteristic time scales of the logarithmic 
susceptibility x(^)- Function x(t) may or may not oscil- 
late in time, but generally x{t) is relatively large within 
a time interval of the order of the relaxation time of the 
system tr [13) 113 • To reveal some qualitative features of 
the effect of vaccination and in particular, its dependence 
on the vaccination period, we will consider s'^^ for two 
types of constraint on this period. 



2. Limited vaccine accumulation 

A different situation occurs if the total amount of vac- 
cine that can be accumulated is limited. This limitation 
implies that ST < M (note that M is the limit on the 
ensemble-averaged amount of the accumulated vaccine). 
Such a constraint is typical for live vaccines, as it may 
be dangerous to store too much vaccine in this case. The 
actual average vaccination rate in this case is now T- 
dependent. We use the notation S^ for this rate, with 
Sa = min(S, M/T). This is S^ that should be used now 
in Eqs. and ([21]) for s^^^ in the limits T > t^ and 
T <^tr, respectively. 
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The behavior of s^^t with varying vaccination period 
T depends on the form of the logarithmic susceptibil- 
ity x{t)- Let us first consider the case where x{i) tias 
a single local minimum (at t = t^), and mono- 
tonically decays to zero with increasing \t — <*|. Here, 
once the vaccine accumulation has reached saturation 
with increasing T (which happens for ST = M), func- 
M\ iBhiXT{t)\ monotonically decreases with 



tion \s 

further increase in T. Indeed, 



j2 j^Mi 



a»T + nT) 



dt 



- a^T + nT) 
(n-a,) > 0,(23) 



V 



^o(t) 



V 



pI/N 



^ R 



t=t. 



FIG. 1: The SVIR epidemic model with susceptible, vacci- 
nated, infected and recovered sub-populations. The arrows 
indicate processes leading to changes of sub-population sizes; 
the corresponding rates are scaled per individual. 



where a* gives the position of the minimum of XT(t) 
over t and is given by equation dxrit* — ai,T)/dat: = 0; 
we choose < a* < 1 and take into account that 
mino<t<T XT(i) < 0. In Eq. (^5)) we have used that, 
if x(i) is minimal for t = t^, then dx/dt > for t > 
and dx/dt < foT t < t*. It follows from Eq. that, 
once the vaccine accumulation has reached saturation, 
further increase in T will only reduce the effect of the 
vaccine. This result is understandable because, if T in- 
creases beyond MfE., the actual average vaccination rate 
Sa decreases. 

A counterintuitive situation may occur if %(<) is os- 
cillating. Here the inequality may be violated. As 
a result, the dependence of the effect of the vaccine on 
T and, consequently, on the actual vaccination rate S^, 
may be nonmonotonic. An example of this behavior is 
discussed in the next section. 



V. RESONANCES IN THE STOCHASTIC SVIR 
MODEL 

We now apply some of our results to an important and 
widely used stochastic epidemic model, the Susceptibles- 
Vaccinated-Infected-Recovered (SVIR) model. The 
model is sketched in Fig.[T] In the absence of vaccination, 
^o(i) = 0, the SVIR model reduces to the stochastic SIR 
model with population turnover, which was originally in- 
troduced to describe the spread of measles, mumps, and 
rubella, see [H, In the SIR model, susceptible in- 
dividuals are brought in, individuals in all population 
groups leave (for example, die), a susceptible individual 
can become infected upon contacting an infected indi- 
vidual, and an infected individual can recover. If we set 
Xi = 5", X2 = /, and — R, the rates of the cor- 
responding processes are: (i) influx of the susceptibles, 
iy(X, r) — iiN for ri — l,ri^i = 0, (ii) leaving, with 
the same rate for all populations, iy(X,r) = iiXi for 
n = -l,rj^i = 0, (iii) infection, W{X,r) = PX1X2/N 
for ri — — l,r2 — l,ri^i_2 = 0, and (iv) recovery of the 
infected, M^(X, r) = 7X2 for r2 = -l,r3 = l,ri^2,3 = 0, 
Fig. ©. 



For /3 > r = 7 -I- /i the SIR model possesses a sin- 
gle endemic state. This state corresponds, in the mean- 
field theory, to an attracting fixed point on the two- 
dimensional phase plane of susceptibles and infected. At 
^ < 4(/3 — r)(r//3)^ this attracting point is a focus. 
The populations of susceptibles, infected and recovered 
exhibit decaying oscillations in time as the system ap- 
proaches the endemic state. It was found in Ref. 
that, in this parameter range, the populations oscillate 
also on the optimal disease extinction path. These oscil- 
lations are illustrated in Fig. [5J 

We will now incorporate vaccination and introduce a 
sub-population of vaccinated X^ = V. The vaccination 
is described by the transition rate M^(X, r) — ^o(i)a^i for 
ri — —1,^4 — l,ri^i_4 = 0. The corresponding term in 
the Hamiltonian Eq. ()12|) has the form ^o(i)-ff'^-' with 

ij(i)(x,p)=xi(ef''-Pi -1). (24) 

Vaccinated individuals leave at the same rate fi as in- 
dividuals in other populations. For simplicity, we as- 
sume that the immunity from the vaccination is never 
lost. In this case fluctuations of the vaccinated pop- 
ulation do not affect fluctuations of other populations, 
and p4 = along the optimal extinction path. Then 
from Eq. (ITil) . the logarithmic susceptibility is x{t) = 

xfjptit)^!! (1 - exp[-p£^pt(t)]) , where xfjp^{t), p^^Jpt{t) 
and xiA are calculated for the SIR model. 

The Fourier spectrum of the logarithmic susceptibility 
x{io) is plotted in Fig.[3](a). It corresponds to the optimal 
extinction path shown in Fig. [5J As one can see, the 
spectrum has a peak at the characteristic frequency of 
oscillations of the system in the absence of vaccination 
ujq (for the chosen parameter values ujq ~ 5.2/i). 

We now consider the effect of the resonant peak in x{uj) 
on vaccination. The dependence of the scaled change of 
the disease extinction barrier s^xt = Q^^^'/N on vaccina- 
tion period T is shown in Fig. [3] (b) . The solid line in 

Fig. [3] (b) shows the behavior of Sg^j where there is no 
limit on vaccine accumulation or, equivalently, for such 
periods where the limitation does not come into play and 
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the actual vaccination rate Sa is independent of T. Func- 
tion |s^xtl = —ScX is seen to be strongly nonmonotonic, 
it displays pronounced maxima (which correspond to the 
minima of Sg^t)- They occur where the vaccination pe- 
riod T coincides with the multiples of the characteristic 
period of the system motion without vaccination 21^/10^. 

For limited vaccine accumulation M, the actual aver- 
age vaccination rate depends on the vaccination period, 
Sa = min(S,Af/r). Beyond a certain value of T, the 
increase of T is accompanied by the decrease Sq. This 
leads to a change of the dependence of Sg^t O'^ ^- 
markably, |s[,^{| still displays resonant peaks at 27rrt/a;o 
with integer n. Their amplitude decreases with increas- 
ing n. The occurrence of the peaks shows that, by tuning 
the vaccination period, the effect of the vaccination can 
be resonantly enhanced; the resonance in this case is in 
the exponent of the disease extinction rate, and there- 
fore it is extremely strong. Counter-intuitively, since the 
actual average vaccination rate decreases with increasing 
T, a strong enhancement of the vaccine can be achieved 
where this rate is decreased. For example, in Fig. |3] (b) 

the maxima of \s^^}t\ fo'^ jiM/'E. = 1 and ^iM/'E. ~ 3 are 
achieved for T in the range where Sq < S. 




FIG. 3: (a) The Fourier transform of the logarithmic suscep- 
tibility in the SIR model. The parameters are the same as in 
Fig. [21 the rescaled frequency is uj' — uj/n- The susceptibility 
spectrum displays a sharp peak at the characteristic vibration 
frequency luq. (b) The change of the scaled extinction barrier 
Scxt ~ M'S^xt/" with vaccination period T. The solid line 
shows s^xt where there is no limit on vaccine accumulation, 
whereas the dashed lines refer to the case of limited accumu- 
lated vaccine amount. The accumulation limit M is scaled 
by the small-T average vaccination rate H, M' = fiMfE, and 
T' — Tfi. The locations of the resonances of s^xt are indepen- 
dent of M. 



0.08 




-^1 

FIG. 2: The most probable trajectories in the stochastic SIR 
model on the plane of the scaled numbers of susceptibles and 
infected, xi = X\/N and X2 = X2/N, respectively. The 
dashed line shows a mean-field trajectory toward the endemic 
state, and the solid line shows the most probable trajec- 
tory followed during the fluctuation-induced disease extinc- 
tion [l||. The plot refers to 13 /fi = 80, and 7/^ = 50. 



VI. CONCLUSIONS 

We have developed a theory of optimal periodic vacci- 
nation against an endemic disease for low average vacci- 
nation rate, as in the case where the vaccine is in short 
supply, or short lived, or cannot be stored in the suf- 
ficient amount. We assume that the vaccination rate is 
insufficient for eliminating the endemic state and thus ex- 
terminating the disease by "brute force" . However, vac- 
cination can change the rate of disease extinction, which 
occurs spontaneously as a result of a comparatively rare 



fluctuation. We show that the optimal vaccination leads 
to an exponentially strong increase of the disease extinc- 
tion rate. This happens because the vaccine changes the 
effective entropic barrier that needs to be overcome for 
spontaneous extinction. 

We find that the optimal vaccination protocol is a peri- 
odic sequence of (5-like pulses. This protocol is essentially 
model-independent, it only requires that the population 
be spatially uniform. In stationary systems, the phase 
of the pulses is irrelevant. In contrast, in periodically 
modulated systems, like in the case of seasonally vary- 
ing infection, it is necessary to appropriately synchronize 
vaccination pulses with the modulation. Moreover, if the 
pulse phase is wrong, vaccination may hamper disease 
extinction. 

For fixed average vaccination rate, the effect of vacci- 
nation in stationary systems increases with the increas- 
ing vaccination period. However, this increase is gener- 
ally nonmonotonic and the disease extinction rate can 
display exponentially strong resonances. They occur if 
the vaccination period coincides with the period of de- 
caying oscillations of the population, which character- 
ize the approach to the endemic state in the mean-field 
(fluctuation-free) approximation. The resonances occur 
also where the vaccination period coincides with a mul- 
tiple of the dynamical period. They are illustrated using 
the well-known SVIR model of population dynamics. 

It turns out that, counterintuitively, the effect of vacci- 
nation can be sometimes enhanced by reducing the aver- 
age vaccination rate. This happens where the mean-field 
dynamics is characterized by decaying oscillations and 
there is a constraint on the amount of vaccine that can 
be stored. In this case lowering the average vaccination 
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rate can allow one to tune the vaccination period in res- 
onance with the system dynamics. 

The analysis is based on the master equation for the 
population dynamics. We solve it in the eikonal approx- 
imation and reduce the problems of the tail of the dis- 
tribution and of the extinction to Hamiltonian dynamics 
of an auxiliary system. A general formulation of the cor- 
responding Hamiltonian problem is obtained for period- 
ically modulated systems. The optimal vaccination pro- 
tocol is found using this formulation with account taken 
of the constraint on the average vaccination rate. The 
feature of the problem that makes it different from other 
problems of optimal control of rare events is that the vac- 



cination rate cannot be negative and it is the average vac- 
cination rate that is given. The analysis can be extended 
to other problems of optimal control of fluctuation-driven 
extinction with similar constraints. 
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